def embed_vid(outvid):
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
return HTML(data='''<video alt="test" width="950" height="500" loop="true" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" loop="true" />
</video>'''.format(encoded.decode('ascii')))
The tracking algorithm is a fork of the Tint (Tint is not Titan) tracking algorithm (http://openradarscience.org/TINT/). The framework has been modified that it can be also applied to model output data that is not not stored in Py-ART radar data container.
#Plot the Medians
um044_n = len(UM044_ens.coords['dim_0'])
um133_n = len(UM133_ens.coords['dim_0'])
cpol_n = len(OBS.dataset['total'].coords['dim_0'])
cpol_n2 = len(CPOL.coords['dim_0'])
CPOL_pdf2 = ravel(OBS.dataset['obs'])
var=['avg_area','dur', 'avg_mean', 'max_mean', 'dist', 'v']
names=['area', 'duration', 'avg-rain', 'max-rain', 'distance', 'speed', '# storms']
#print(CPOL_pdf[var].median())
medians = pd.DataFrame({'a': list(CPOL_pdf[var].median().values)+[int(cpol_n)],
'b': list(CPOL_pdf2[var].median().values)+[int(cpol_n2)],
'd': list(UM044_pdf[var].median().values)+[int(um044_n)],
'c': list(UM133_pdf[var].median().values)+[int(um133_n)]} )
#index=('Area','Duration', 'Mean-Rain', 'Max-Rain'))
medians.index=names
medians.columns=['Bootstrap', 'Cpol', 'UM 1.33km', 'UM 0.44km']
#print('Medians:')
medians.round(2)
| Bootstrap | Cpol | UM 1.33km | UM 0.44km | |
|---|---|---|---|---|
| area | 111.20 | 119.38 | 81.64 | 67.34 |
| duration | 60.00 | 74.67 | 60.00 | 60.00 |
| avg-rain | 4.52 | 4.50 | 4.91 | 5.93 |
| max-rain | 6.58 | 6.60 | 7.37 | 9.06 |
| distance | 15.10 | 16.42 | 10.60 | 12.09 |
| speed | 12.61 | 12.83 | 10.03 | 12.78 |
| # storms | 1184.00 | 42.00 | 50.00 | 73.00 |
#Create Hex-Bin plot
mpld3.disable_notebook()
var=['avg_area','dur', 'max_mean', 'avg_mean']
medians = pd.DataFrame({ 'Bootstrap':CPOL_pdf[var].median(),
'CPOL':CPOL_pdf2[var].median(),
'UM 1.33km':UM133_pdf[var].median(),
'UM 0.44km':UM044_pdf[var].median()})
#medians.loc['avg_area'] /= 2.5**2
histdata = [UM044_pdf[var].dropna(), UM133_pdf[var].dropna(), CPOL_pdf2[var].dropna(), CPOL_pdf[var]][::-1]
titles = ['UM 0.44km ens', 'UM 1.33km ens', 'CPOL', 'Bootstrap'][::-1]
fig = plt.figure(figsize=(10,8))
colm = colmap2.Blues
colm.set_under('w', alpha=0)
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, sharey=True)
fig.subplots_adjust(bottom=0.07, right=0.98, left=0.01, top=0.35, wspace=0.01)
cbar_ax = fig.add_axes([0.15, 0.02, 0.7, 0.01])
hexbin_data = []
gridsize = 10
vmin=0.0001
vmax=0.1
nticks=10
YMax, XMax = 260, 70*2.5**2
#YMax, XMax = 50, 70*2.5**2
for i, ax in enumerate((ax1, ax2, ax3, ax3)):
data = histdata[i][var[:2]]
#data[var[0]] /= 2.5**2
#data = data.loc[(data[var[0]] <= XMax) & (data[var[1]] <=YMax)]
X = data[var[0]].values
Y = data[var[1]].values
ax.set_ylim(0,YMax)
ax.set_xlim(0,XMax)
hb = ax.hexbin(X, Y, gridsize=gridsize, extent=(0,XMax,0,YMax))
hexbin_data.append(np.ones_like(Y, dtype=np.float) / hb.get_array().sum())
cl = plt.cla()
medians = OrderedDict()
for i, ax in enumerate((ax1, ax2, ax3, ax4)):
data = histdata[i][var[:2]]
#data[var[0]] /= 2.5**2
#data = data.loc[(data[var[0]] <= XMax) & (data[var[1]] <=YMax)]
X = data[var[0]].values
Y = data[var[1]].values
ax.set_ylim(0,YMax)
ax.set_xlim(0,XMax)
im = ax.hexbin(X, Y, gridsize=gridsize, cmap=colm, marginals=False, extent=(0,XMax,0,YMax),
vmin=vmin, vmax=vmax, C=hexbin_data[i], reduce_C_function=np.sum)
ax.set_title(titles[i], fontsize=24)
ax.grid(color='w', linestyle='', linewidth=0)
ax.tick_params(labelsize=24)
ax.xaxis.set_ticks(ax.xaxis.get_ticklocs()[:-1])
x, y = histdata[i][var[0]].median(), histdata[i][var[1]].median()
z, zz= histdata[i][var[2]].median(), histdata[i][var[-1]].median()
sx, sy = histdata[i][var[0]].std(), histdata[i][var[1]].std()
medians[titles[i]] = '%2.1f km^2, %2i min (max: %2.1f mm/h, mean: %2.1f mm/h);'%(x,y, z, zz)
ax.hlines(y*np.ones_like(histdata[i][var[0]]),0,x, 'firebrick', lw=3)
ax.vlines(x*np.ones_like(histdata[i][var[1]]),0,y, 'firebrick', lw=3)
ax.scatter([x], [y], marker='o', s=400, c='firebrick', alpha=0.8)
if i == 0:
ax.set_ylabel('Duration [min]', fontsize=24)
ax.set_xlabel('Rainfall Area [km$^2$]', fontsize=24)
ary=im.get_array()/im.get_array().sum() * 100
im.set_array = ary
cbar = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
tmp = cbar.ax.tick_params(labelsize=24)
tmp = cbar.set_label('Density [ ]',size=24)
#cbar.set_ticks(np.linspace(vmin, vmax, nticks).round(2))
#cbar.set_ticklabels(np.linspace(vmin, vmax, nticks).round(2))
#print('Medians:')
#medians
#print(' '.join(['%s: %s'%(k, v) for (k,v) in medians.items()]))
#plt.show()
<matplotlib.figure.Figure at 0x7f69f09e0e80>
colm = matplotlib_to_plotly(CubeYF_20.get_mpl_colormap(N=8, gamma=2.0),8, rgb=False)
#Plot the tracks
mpld3.enable_notebook()
fig = plt.figure(figsize=(14,5.5))
fig.subplots_adjust(right=0.94, bottom=0.45, top=0.95,left=0.01, hspace=0, wspace=0)
#cbar_ax = fig.add_axes([0.12, 0.17, 0.74, 0.02])
o = namedtuple('Sim', 'dataset percentiles')
tmp_obs = o({'obs': CPOL_t.dataset['obs']}, {'obs': CPOL_t.percentiles['obs']})
with nc(CPOLF) as fnc:
lon=fnc.variables['longitude'][:]
lat=fnc.variables['latitude'][:]
tp = 'mean'
num=80
titels = ['CPOL', 'UM 1.33km', 'UM 0.44km']
m = None
for i, data in enumerate(((tmp_obs, storm_obs), (UM133_t, storm_UM133), (UM044_t, storm_UM044))):
tracks, stormId = data
ax = fig.add_subplot(1,3,i+1)
#ax2 = fig.add_subplot(2,3,i+4)
ax.set_title(titels[i], fontsize=18)
for ii, tr in enumerate(tracks.dataset.keys()):
if ii == 0:
draw_map = None
m2 = None
else:
draw_map = m
Id = [stormId.Sim[tr]]
perc = tracks.percentiles[tr][tp][num]
ax, m, im = plot_traj(tracks.dataset[tr], lon, lat, ax=ax, mintrace=2, thresh=('mean', perc),
color=colm[ii], create_map=draw_map, basemap_res='f', lw=0.5, size=20, particles=None)
#ax2, m2, im = plot_traj(tracks.dataset[tr], lon, lat, ax=ax2, mintrace=2, thresh=('mean', perc),
# color=colm[ii], create_map=m2, basemap_res='f', lw=0.5, size=20, particles=Id)
#break
#plt.show()
mpld3.disable_notebook()
variables=dict(dur=('Duration [min]',300), avg_area=('Area [km$^2$]', 100*2.5**2), avg_mean=('Rain-rate [mm/h]',12))
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.6,left=0.01, hspace=0, wspace=0)
i = 0
for y, prop in variables.items():
label, ylim = prop
npl = len(list(variables.keys()))
i += 1
ax = fig.add_subplot(npl,1,i)
ax = sns.boxplot(x="quant", y=y, hue="run", data=df_stack, palette="muted", ax=ax)
#ax = sns.stripplot(x="quant", y=y, hue="run", data=df_stack, jitter=True, palette="Set2", dodge=True)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
ax.set_xlim(0.5,5.5)
ax.set_ylim(0,ylim)
ax.yaxis.set_ticks(ax.yaxis.get_ticklocs()[1:])
ax.set_xlabel('Quintiles', fontsize=26)
ax.set_ylabel(label, fontsize=26)
#break
#fig.savefig(os.path.join(os.environ['HOME'], 'Todds_Rainfall_2.pdf'), bbox_inches='tight', format='pdf', dpi=300)
P = list(np.arange(0,99))+[99, 99.9, 99.99, 99.999, 100]
PERC = pd.DataFrame({'Obs':np.percentile(CPOL_pdf2['avg_mean'].values, P),
'UM 1.33km': np.percentile(UM133_pdf['avg_mean'].values, P),
'UM 0.44km': np.percentile(UM044_pdf['avg_mean'].values, P)},index=P)
from mpl_toolkits.axes_grid.inset_locator import inset_axes
#Plot the percentages
fig=plt.figure()
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.plot(PERC.index, PERC['Obs'].values, color='lightseagreen', linestyle='-', label='CPOL',lw=3)
ax.plot(PERC.index, PERC['UM 1.33km'].values, color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
ax.fill_between(PERC.index, member.min(axis=0)-0.25, member.max(axis=0)+0.25, color='cornflowerblue', alpha=0.2)
ax.plot(PERC.index, member.mean(axis=0), color='lightseagreen', linestyle='--', label='Bootstrap',lw=3)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile [\%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
ax.grid(color='r', linestyle='-', linewidth=0)
#Bootstrapping stuff:
Cold-Pool tracking using density potential temp. pertubation fields
embed_vid(os.path.join(dest_dir, 'ColdPool-Ens-2.mp4'))
mpld3.disable_notebook()
variables=dict(cp_area=('Cold-Pool Area [km$^2$]', 220), cp_mean=('Cold-Pool strength [K]', 2),
cp_dur=('Cold-Pool Lifetime [min]', 120))
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.6,left=0.01, hspace=0, wspace=0)
i = 0
for y, prop in variables.items():
label, ylim = prop
npl = len(list(variables.keys()))
i += 1
ax = fig.add_subplot(npl,1,i)
ax = sns.boxplot(x="quant", y=y, hue="run", data=df_stack, palette="muted", ax=ax)
#ax = sns.stripplot(x="quant", y=y, hue="run", data=df_stack, jitter=True, palette="Set2", dodge=True)
if i == 1:
ax.legend(loc=0, fontsize=24)
else:
ax.legend(loc=0, fontsize=0.001)
ax.tick_params(labelsize=24)
ax.set_xlim(0.5,5.5)
if i == 3:
ax.set_ylim(60,ylim)
else:
ax.set_ylim(0, ylim)
ax.yaxis.set_ticks(ax.yaxis.get_ticklocs()[1:])
ax.set_xlabel('Quintiles', fontsize=26)
ax.set_ylabel(label, fontsize=26)
#break
#fig.savefig(os.path.join(os.environ['HOME'], 'Todds_Rainfall_2.pdf'), bbox_inches='tight', format='pdf', dpi=300)
mpld3.disable_notebook()
variables=dict(omega=('Vert. motion', (0, 6)), bflux=('Bouyancy Flux', (-3,30)), mfluxdiv=('Moist. Flx Conv',(-5,6)),
momflux=('Momentum Flux', (-2, 2)), moistflux=('Moisture Flux', (-0.01, 0.04)))
#cloud_w=('Cloud Water', (0,0.005)))
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.6,left=0.01, hspace=0, wspace=0)
i = 0
for y, prop in variables.items():
label, ylim = prop
npl = len(list(variables.keys()))
i += 1
ax = fig.add_subplot(npl,1,i)
ax = sns.boxplot(x="quant", y=y, hue="run", data=storm_prop, palette="muted", ax=ax)
#ax = sns.stripplot(x="quant", y=y, hue="run", data=df_stack, jitter=True, palette="Set2", dodge=True)
if i == 1:
ax.legend(loc=2, fontsize=24)
else:
ax.legend(fontsize=1)
ax.tick_params(labelsize=24)
ax.set_xlim(0.5,5.5)
ax.set_ylim(*ylim)
ax.yaxis.set_ticks(ax.yaxis.get_ticklocs()[1:])
ax.set_xlabel('Quintiles', fontsize=26)
ax.set_ylabel(label, fontsize=26)
fig = plt.figure()
from matplotlib.ticker import Formatter
from scipy import interpolate
variables=dict( momflux=('Momentum Flux', (-1.5, 1.3)),
mfluxdiv =('Moist. Flx Conv', (-0.5e-3,1.3e-3)),
mflux=('Moisture Flux', (-1.2e-4, 1.6e-3)), omega=(('Vert. Motion'), (-9.5e-2, 8.5e-1)))
i = 0
fontsize=24
y = np.linspace(P[-1], P[0], 40)[::-1]
for x, prop in variables.items():
label, xlim = prop
npl = len(list(variables.keys()))
for ii in range(1, 6):
i += 1
ax = fig.add_subplot(npl,5,i)
for ff, fluxes in enumerate((fluxes133, fluxes044)):
flux = fluxes[x][ii]
f = interpolate.interp1d(P, flux, kind='slinear')
#xx = f(y)
y = P
xx = flux
ax.plot(xx, y, label=('UM 1.33km', 'UM 0.44km')[ff])
ax.set_ylim(max(y), min(y))
ax.set_xlim(*xlim)
if ii == 1:
ax.set_ylabel('Pressure [hPa]', fontsize=fontsize-2)
else:
ax.yaxis.set_ticklabels([])
if i == 1:
ax.legend(loc=1, fontsize=fontsize)
else:
ax.legend(loc=1, fontsize=0.001)
ax.set_xlabel(label, fontsize=fontsize-2)
ax.tick_params(labelsize=fontsize-2)
ax.ticklabel_format(axis='x', style='sci')
if i < 6:
ax.set_title('Qantile %i'%ii, fontsize=fontsize)
#ax.xaxis.set_major_formatter('%.1E')
ax.ticklabel_format(style='sci', axis='x')
ax.xaxis.major.formatter.set_powerlimits((0,0))
fig.subplots_adjust(right=0.98, bottom=0.25, top=0.99, left=0.01, hspace=0.15, wspace=0.04)